{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 360x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEKCAYAAABT352BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hU1dbA4d+aFJIQeg0lhKJIUVqQXkIRpIoUEbyickVABJFProBKpCggiMLFCNgF0YgoIBcVIYBIEVSaIL0FAoReQs/6/jiTkEAqM5OZhP0+zzwzc+bsc9YEXTnZZ++1RVUxDMMwXMfm7gAMwzByOpNoDcMwXMwkWsMwDBczidYwDMPFTKI1DMNwMW93B+BMhQsX1pCQEHeHYRhGDvPHH3+cUNUid9o+RyXakJAQNmzY4O4wDMPIYUTkgCPtTdeBYRiGi5lEaxiG4WIm0RqGYbiYSbSGYRguZhKtYRiGi2VpohWRj0XkuIhsTbKtoIgsEZFd9ucC9u0iIlNEZLeIbBaRmukdf8sfWyhd+igjZkYx4bcJae474bcJjJgZRalSMYRIiGnnYDuA2bNJ1m727HSbJGsn0sS0c2M7gJiYGJo0acLRo0cz3siBdncNVc2yB9AYqAlsTbJtAvCK/fUrwHj76zbAYkCAusC6DBxfCeygDC2sw2cs07QMn7FMGVpYCeygNtPO4Xbdp47XXPctU+hntaOf5rpvmXafOj7D7cCW4XazZqkGBGiydgEB1nbTLvPtEtoGBlrtAgP7ZahNgn79+qnNZtN+/fplvFE2AmxQB3KfaBaXSRSREOAHVa1qf78DaKqqMSISBCxX1YoiMt3+es6t+6Vx7MQv4y3CM/fey0VvX2Zt3QxA7/pNiI89QpOgIJ5dtYrrKXx3G8JHjRoCsOzIEfwKF2fG2lUA9HygBl9t2UQ8t7cT8aNrxfIUkXhCCxcGYHF0NMUqVGbKzz9gs/mjejmF89m4oTcAeKRsWe4JDKRKgQIAfLd/Pz8cOkI8N25r5y1Cm1KlqNmsDSM//YCDew7yQlhDHihYkPJ586b6/bwQPmzUkPPXrvHL4cM0e+wZBr0dzq+LlzPxuSdZdCiaGyl9P3zpXPoemhe7l1L+RYi9fJmlR3fy1PDXaNG3DTabD6rXb2uH+KDxV5k//28mTfqa0NDnCQwsxrFjm9m+fS6/rhoHei2Fdt40ajiMjz/+PypUyMuHH67l88//R4MGr+DjE8Do0bmAq7e3w4/XXrvE/v1R7N8fxY8/jiIgAMaNW8LixatZuXIccPu/A/gQFvYGDRsOA2DHjvmcPr2Tn39+2QpHfIHb4xTxIz7+EgMHfsWuXSeoXXsAAFu2zMbL6xzz5r2U4r87+NCjxwfMnv0MAE8+OZOzZ32pVq0XAKNH+wAp/Dzt32/duveoUKEE06Z1BaBjx3fw9Q1h7tyeqXw/b8aM+Z6WLduSJw98/PEYGjQI5ZFHWgPg5eVDfPzt5/Px8SMu7jyjRo2iRYsWNG7cmCtXrjB27FhatWpFixYtuHz59vP5+flx6dKlFOLInkTkD1UNveMDOJKl7+QBhJD8ivbMLZ+ftj//ADRMsn0pEJrC8foAG+wPFdCeoEdA24FWtZH4WynUz0frg94APQxaFOxXX2gAaGHQJvbPb4DWAm2Yyyux/T1eom1Ae9j3B9TLfr4hDNFSgvZK0r4EaOu8AaqqOoSX1BfUO8n5vEG7+PgmHt8XdFiS9gLaw9dfe4D629v52M+30/75YyFlVFV1/cp1KqDv2dv+Yd/f1/7sZ3+ebP98u739s7Vqq6rq9NGTVEBn2L9fLvv+ueznG5CrlgroYnv7pfb2r3Topqqq/XyrJu6f9HmIb6iqqg5s000F9MPC/jqrQkF9vHReBTRfvqlKQHPFZu2PN0pQRaXAswqiS5ZEq6pq585TFEThhIqoUqq5tT/+9mcf67nWK9bnhCuInjhh/Wzr1x+u4KtwRKGHgre9XYBCT4VnFAqpiNrb91EokeSKppt93wB7Oy+FPAoxqqpasmRnhfuTtG+vvr61kpzPdsv56mnu3E0Tj58nT32Fh5K0f0ChZJLziUJphRj75/dq8eI9Ett7e4coPJXkfNxyvgIKLyio/ZFHK1cekuT7+ShUSXI+FO7XUqVi9MqVKyoiOnbsWFVVPXfunIqITpw4UY8cOaKdO3dObBMQEKCPP/647tu3T3MSHLyi9eREuyiFRFsrrWNLwn8ggR20TJm0f3BlyqgS2D5ZEjLtUm9Xuvw5XX1wtf60+yed+/dc/eSvT3TK2im659QeVVUNqvGnEnJvsnZBJYvqfVVXqqrqX5FTdEWdIN1asYAeKRqgl31sqqCN7l2s0DexnQ3034LG5EZ3FkT3bVujqqp/vvuKbm5WVXd3bqbHnntCJ5TqofcHeSuI/XyiVM+lxeqk3cVRpozaz2dT8Ev88zpDP5ds3C4oSHXFCtUfflCdM0d1xgzVpUutNtevJyTf29uJqF69qnruXOrn7Nu3r9psNvXz81ObzaYtWrTQwoUL686dO9MONhvJCYl2BxBkfx0E7LC/ng48ntJ+qT0qgxLYUbnHN2N9kff4KoEddaNp53C7WbNUbYENleq5rHbVc6ktsGHq/Xzx8aqnTunmR17TTqD9QTfanzuBRgdX0J0P1dLrZ06rquqCfs11Z0ErAV/wsS7LUmq3qdPI9OO0dVLor7BRob/abJ0y1PeZk9tZCTp5O+ikZcqo/vyzqq+vauvWqhERqtHRyduGhnbSwECrXWBgf61Uqan27dtXb9y4oaqq33zzjc6bN0/j4+PTDsKD5YRE+zbJb4ZNsL9uS/KbYb+nd+wggrRMGStZjF+V9s2U8avG6/AZy7RMGdVwRpp2DrZbtneZ5hlVWIvVWabhjNRidaz3y/amnaCTtlNIs13c1TjddHSTRm6N1DHL3tBOb/bVfI+8pgp6JLeXnvETnft0Xf1x63y9dO1SquecNctKLCLWc0Zv+uTkdjdvot18JNxE++cf1ZdeUi1f/uZntWpZCTetdgkaN26sjRs3Tnx//vx5VVU9cuSINm7cWGNiYjL2xdwoWyVaYA4Qg3VXIRroDRSydwvssj8XtO8rwDRgD7Alpf7ZWx+1atVy8o/XyKjxq8bflhyX7c1YYk9sBxlul7Bf4QmFVUFDBwXo6gcKqILuyY8+0T2XDl8yzGlx3g3SS9Dx8arbtqmOG6fapo3V5WBdCd/+SNrFce3aNT18+LCqWkm2QIECOnny5Gw1UiFbJVpXP0yizeZGjszwrglJdtneZaojRya+X//Jm3q2Yogq6LcfWjd7zl85r80+a6YTf5uon2/8/Ga7W49jZJp1Y+72h0jK+586dUq9vLyS3HC7+fDz88va4DPBJFqTaO9KaV6ZXr9u3fWxi54+UVuPqaSEo4Sjwe8Eq/8Yfx3wvwEmyTooI1e0tzpy5Ij26NFDAwKsEQ42m0179Ojh0V0IjiZaMwXXyJaGNhhKWNmwZNvCyoYxtMFQ8PKCtm2tjRcuUPI/Y1g8ei+nz/QjotF47i18L1dvXKXQ+P/SL7TfbccxMm7sWAgISL7N29vanpqgoCDy5s3L5cuX8fHxQVXJmzcvxYsXd22wbmQSrZGzBQbC5s3QrRv5342gb/eJ/PdAVQr55CN8BURsiKD/ov7sPLnT3ZFmSz17wowZUKYMiEDevNC/v7U9LceOHaNv376sX7+efv36cezYsawJ2E2yfGaYK4WGhqpZYcFI1YYNnHm+N4EbNvPHT59Sp+VTfLdtHp0jO+Nt8+at5m/xYt0X8bJ5uTvSHCE2FopkcPGXWbNm8ddffzFp0iTXBnWHHJ0ZZq5ojbtHaCh/VymCdzzUafkUAJ0qP0p8uDJ5TT7+b8n/0fCThvxz4h/3xpkD/PgjhITA0qUZ23/btm2sW7cuxem8OYG5ojXuXiLWc7t26GefMefwj7yw+AVsYuPAiwcI8AlIu72RqjNnoHFj2L8fVq6E6tXT3v/q1at4e3tjs3nmtZ+5ojUMR0ydCj/9hNSqRY8r9/J3/7+Z1WkWAT4BqCr7Tu9zd4TZUv78sHix9fzww1bCTYuvry82m43z58+zaNGiLIkxK5lEa9y9Ro6EAQPg118hPh4aNKD4Vz/QqkIrAL7Y/AX3TbuPt359i+spVLYy0laypJVsL1+G1q3h3Ln027z++us8+uijOa6urUm0xt0rPNx6rlMH/vwTmjcHX9/Ej1tXaE3Hih0Zvmw4Ie+G8NGfHyVrHrUvYwXR72ZVqsCCBdC9O+TJk/7+I0aMYPny5TluqJfpozWMBKo3+23nzYP77oPKlZm7bS7/XvBvzl45ywsPvsCUh6cQtS+KbnO7Edkl0ozDzYRDh6BECWuoc3ouXrxI7ty5XR9UBpg+WsNwloQke/UqDBkCtWvDl1/SpXIXdg/cTVhIGJ9t+ozXo143SfYOHD8ONWvCCy9Yv9PSMmfOHEJCQoiOjs6a4FzMJFrDuJWvL/z2m5UVevaE/v0p7JWHZb2WMajOIEavHE2dknXI75ff3ZFmK0WLQu/eEBEBb72V9r516tShTZs2+CbpysnOTNeBYaTm2jUYMQLefhtCQ1nxxRi6LHyCOVsq0qrCGrzFm887fc5jVR9zd6TZRnw89OoFs2ZBoUJw6hQEB1tTdtObTeZOpuvAMFzFxwcmTIDvv2dfvcp0WfgEkV0iafHFb0R2iSSeeLp/250RS0cQr/HujjZbsNmgRQvr+eRJqwvhwAHo04cUV+uNjo6md+/enD17NuuDdSKTaA0jPR078s1jVaw+2WP+AHSu0J4fHv+B2iVq8+aqN+kc2Zmc9NehK40caV3ZJhUXZ/3xcKuYmBgiIyNZv3591gTnIt7uDgBARAYBz2IV+56pqu+KSEHga6wVGfYD3VT1tNuCNO5qQ5fEwRvNbm7IlYtWwEOvv860tk9y7cY1JOFmmpGmgwczvr127docOnSI/Pmzd3+4269oRaQqVpJ9EKgGtBORe7CWtVmqqvdgrbzwivuiNO564eE3y62C1cHo74+UKMGA2s8zuN5gAH7e8zNL9ixxX5zZQHBw5rYnJNlVq1Zx9WpKS8x7PrcnWqASsFZV41T1OrAC6AR0BD6z7/MZ8Iib4jOM223eDA0bQt++MHcuYBXRH/vrWFrPbs27a981XQmpSKmGLVj9tKnZtGkTjRo1IiIiwnWBuZAnJNqtQGMRKSQiAUAboDRQTFVjAOzPRd0Yo2HcNHKkNer+xx/hyy/h0UcBkLNnWdRjER0rdmTwT4N5ZsEzXLl+xc3Bep5ba9iWKmXVsf36a2sIc0qqVavGrFmzePbZZ7M2WGdxZHkGZz2wFmn8E1gJfABMBs7css/pVNr2ATYAG4KDg+9snQrDcFRMjGrRoqpDhuiNS3EaHhVuLZszOVh/2PFDsl3NYpC3W7jQ6pcZMSL9fQ8ePJjlq+eSE5ayUdWPVLWmqjYGTmGtiHtMRIIA7M/HU2k7Q1VDVTW0SEarDBuGs+XNC507w6RJ2OrVZ2Sxbnzb7VuqFavGU/OfImpfFEDi1N3aJWq7OWDP0q4dPPUU/PILXE+jfk9sbCw1a9bk119/ZdSoUVkWn6M8YsKCiBRV1eMiEgz8DNQDhgMnVXWciLyCtQz50LSOYyYsGG63cCE88wxcvAjvvAPPPUfU/uV0juzMgyUf5I+YP8zU3VRcvGhNyvPxSflzf3//FAuD+/n5cenSJZfGllMmLHwrItuAhcDzag3jGge0FJFdQEv7e8PwbO3bWzfKGjWCJdbog7CyYVQoWIG6H/6Et3ibguKpyJ3bSrKnTsFXX93++d69e+nRowcB9jtpAQEB9OzZk337PL9msEckWlVtpKqVVbWaqi61bzupqs1V9R778yl3x2kYGRIUZBVi/eILEGHtsi8I/n0H4SvgeNxx6n5Ul74/9OXUJfOfdErGjYMePaxyE0klXT034So2u6ye6xGJ1jByHJsNAgKI2hfFwZefZe6HVtXrX+p9gL+3PzP/mEnE+uw5VMnVXnvNWm+sVy+rOyGphNVzBw8ejKqyc2f2WL3YI/poncX00Rqe5rdnWtLgk19u276iY3Ue/GY1/j7+rDq4ivx++alatKobIvRMK1ZA06bWAhhTp97++alTp/j9999p2bIlXhkpbuugnNJHaxg5UoOPlySfUTZwIPj60qRmJ/x9/FFVBv80mBrTazB0yVAuXL3g3oA9RJMmMGgQ/Pe/Ka+kW7BgQVq3bp0lSdYZTKI1jKz03nuwZw+8+CIAMn8+K/+sxsDSXXl79dtUnlaZ77Z/Z2aVAW++aU28q1Qp5c8vXbrExIkTiYqKytrA7oBJtIaRVUaOtJ4TpkIB/P03/jM/ZVK/74iOfoxy1wJ5NPJRnl3wbOLY2wR32xplAQFWkfASJVJekcHHx4dJkybx448/Zn1wmWQSrWFklYTFIJMaMQJ27IDHHqPkx98Q9cYB1l98gsfvf5xuc7vxwYYPuPH6a3f1RIfjx611MxcvTr7d29ubbdu2MX78ePcElgnmZphheIqdO+GNN6BVK3jySRZtmsu/5nTl1HgoMqHwXTvR4coVqFULTp+GrVuhQIHb91FVl5apNDfDDCOnuPdea5mBJ58EoO3SgxycapWMLpe/HE1DmroxOPfJlQs++wyOHbNukN1q8uTJhIWFeXS/tkm0huGJwsNhyBAC46yJ/+v6/I7YbCl3P9wFatWyelm++MJa5NFms8bazp4N+fLlo0SJEi6fhusIk2gNwwNF9WpCkQmFWbvsCwC2lcuD76sw6aE8bo7MfcqVs8oqxsYmX2ssV65n+PLLLxOn5noik2gNwwOtP7KeyC6R1A17AoDKe8+zeMsDXIu/5ubI3GfkyNtHHyRda+zIkSNcueKZ9X/NzTDD8HTh4dZqho0aQcuWHL1wlBNxJ+66mWQ2W8rDvETg9983UKdOHb7++mu6dOni9HM7ejPMIxZnNAwjDbf0y/b8qhvbzu5iTe81hOQPcUtI7hAcbHUXpLS9Ro0ajBo1itDQO86FLmW6DgwjO5k0iYXvnyH+8iVaz2rNybiT7o4oy6S01pi/v7Xdy8uLESNGEBIS4pbY0mMSrWFkJ+XLE/DXFv7a1Zz9Z/bTfk57Ll3z3LvtznTrWmNgVfjq2fPmPuvWrePnn392T4BpyFSiFZFcInKPiISKSDlXBWUYRioeeQQGD6bEZ/NYnvt51kavZdjSYe6OKsv07An798ONG9YohFurJL700ksMG+Z5P48MJVoRaSYii4DTwD/A78AuETklIh+KyD2OBCEig0XkbxHZKiJzRMRPRMqKyDoR2SUiX4uIryPnMIwcY9w4qFOHuiNnsujB93it8WvujijLiVhrjC1bZiXeBB9//LFHFplJM9GKSICIfA48B8zHWgq8KlABeAB4DNgOTBORFOZspE9ESgIDgVBVrQp4Ad2B8cBkVb0HK8H3vpPjG0aO4+trrc0dEMDDpwpSKKAQV65f4ft/vnd3ZFnq3/+2FnMMDr65rWLFiuRNKNjjQVJNtGJNHA4HJqjqY/bVZper6nZV3auqf6vqElWdpKoPAZtFpGdqx0uHN+AvIt5AABADNAPm2j//DHjkDo9tGDlPmTKwe3diB+V7696j09ed+OSvT9wcWNYJCrKKzdhuyWJr1qyhU6dOHjVTLK0r2lBgvKpuzciBVDUK2CYigZkJQFUPAxOBg1gJ9izwB3BGVRMWHo4GSqbUXkT6iMgGEdkQGxubmVMbRvYWaP9fbf58Bh8OpmW5lvRe0Jtxq5KvY5qTyyueOgUvvQSrV9/cdunSJf7880/27NnjvsBukeo4WlVdn9J2EckNNAH8gI2qujdJm78yG4CIFAA6AmWBM8A3wMMphZRKnDOAGWBNWMjs+Q0jW4uPh4kT8dm4kXmrV1Dj9GMMWzqMfLny0a92v8TyipFdIt0dqUv4+cFHH8HJk1C/vrWtadOm7Nu3D9utl7pulNlRBw8Am4A3gCHAchEZ6mAMLYB9qhqrqteAeUB9IL+9KwGgFHDEwfMYRs5js8GcOeDnR+ATT/Pr4z9TLHcxXlj8Aq8uezUxyebU8ooBAfDYYzB3Lpw/b22z2WzYbDZUlatXr7o3QLv0boY9dMumZkAlVa2tqg1UNRhwtMrFQaCu/cabAM2BbUAUkDCXrhfWzTjDMG5VqpRV1mrzZooPf5OVT6/k3zX/zdhfx9IvtF+OTbIJnn7aqnnwzTc3t128eJFKlSrx9ttvuy+wJNK7oi0uIm+IiI/9/UXgXRHpJyK9ReQNrKvNO6aq67Buev0JbLHHNAP4D/CSiOwGCgEfOXIew8jRWreG4cPhww+58Mtivt3+Lcv2Nub99e/ftiROTlO3LlSsCJ8kuQ+YO3duHn74YapUqeK+wJJIt6iMiFTASnqTVPUfEekBPIo1OmAt1hCs8y6PNANMURnjrnb9Ols+HkezM+8S2fUbwso144H37yfmQkyO7j4AmDIF/vwTPvwQvF1QwcXlRWVUdbeI9ANGiMgRVZ0JfHmnJzQMw0W8vVlcxZfIEt8QdiUIgC3HtzAmbAzrj6zP0Yl24MCUt8fFxbF582bq1q2btQHdIkM3w1T1uqq+AewVkWkiUtDFcRmGcQeGLokjrFyzxDW6NRxGNHmV//v5gnsDyyKbNlkDMRK89NJLtGzZkosXL7ovKDKQaEWkiIjUFJEAVV0KvA6MFZFmrg/PMIxMCQ+3irbaV4b9e2APJBw+7hicZrOcYNEiqF7dmpabYODAgfzwww/4+/u7LzDS6aMVkf7Am1g3wRRonTCBQUT6AsWBMUkmFriV6aM1DDvVxClT/3nhPn6u4sefff506Uqx7nb5sjVbrE0bay0xZ3L1KrghQGFVLYlV2yCxk0dVP8CaXDDzTk9uGIaLiFijEEJDeeuTQ/za4KMcnWTBmrzQowfMmwdnztzcfuLECUaPHs3+pNVnslh6ibYA8ICIlMWakpsr6Yeq+jfQ10WxGYbhiLFj4fvvsT0/gMCK93M9/joXrubsvtqnn7aubCOTTISLi4sjPDzcrVW90us6qI91xVoOWAb8S1VPZVFsmWa6DgwjZVdvXCVscjVCKzbjvfbT3B2Oy6jC/fdDsWKwdOnN7ceOHaNYsWJ3fFyXdR2IiJ+qrlbVKqrqr6ptM5Jkk0ybNQzDQ/heuMT8tw9T4c0Idp/a7e5wXEbEuppdsCD59mLFihETE0OTJk04evRolseVVtfBffYr2gwTkV4OxmMYhivky4d/18d5Ya2yZMTj7o7GpSpXhty5b9/eqlUrVq5cyahRo7I8plQTrapuxEq2L9krdqVKREJE5APgd08ZgWAYRnK5353GntBy9J6xgb+/m+HucFxqwQJo2RKuXwd/f39EhC1btgAQERGBiGTpkK80b4ap6sfAXmCLiCyzT1Z4Q0ReE5HxIvKliGzBqkPwpqpuz4qgDcO4A97eFF0YRXRBL0o/PQgOHXJ3RC5z44a1+sKSJbB371569OhBgH0J3YCAAHr27Mm+ffuyLJ50Jyyo6vdAZaxVDgpjrXTwBNAAOA70V9XmqnrQlYEahuG4PMWDkfkLyNOohTUeKodq2xYKF7YKzQQFBZE3b14uX76Mn58fly9fJm/evBQvXjzL4snQjStVvYyVaD9zbTiGYbha2fptYGEbzl4+S+4rl/D29bu5fncO4etrrfITEWGtwnDs2DH69u3L/v372bhxY5bfEPOcEuSGYWSZvaf3cv/EchxvUCNxum5O89RTcPWqVRd93rx5TJs2jUceeYRBgwYxb968LI3FDMUyjLtQ2fxlKVX8XtbFb+KR4cORKlWgfXt3h+VU1avDc89BhQo3tz377LNuicUjrmhFpKKIbEzyOCciL4pIQRFZIiK77M8F3B2rYeQEIsKEh96m58OXOHpPkPV39rZt1ofh4W6NzZk++ABatUq+7dq1axw7dixL47jjRGsf0hXijCBUdYeqVlfV6kAtIA74DngFWKqq9wBL7e8Nw3CChsENeej+jjR75Czx/n7QsSOcPg1vvOHu0Jzq6FFIOvu2fv36PPXUU1kaQ4YTrYj8T0TGiUhR+8SE3cDv9ipeztQc2KOqB7BWx024AfcZ1ogHwzCcZFyLcezKfZn5b/ayqn1l8ZVeVnj+eejeHa5ds94PGTKEfv36ZWkMmbmiva6qrwCXgMnAFFUtCgQ5OabuwBz762KqGgNgfy56684i0kdENojIhtjYWCeHYhg5232F72P789vpdCg37NyZWDAcEeuRA7oRnnoKjh+HxYut9927d6dDhw5ZGkNmEm1COd1BwBXgVft7pw3GExFfoANW+cUMUdUZqhqqqqFFihRxViiGcde4p9A9EB5O7IXjN+sLHjtmVWjJAYm2dWuryEzSxRsPHDjAjh07siyGzCTaIiISgZVgX1bVOBF5GHjOifE8DPypqgl/vxwTkSAA+/NxJ57LMAy7Z+Y/Q4l3SrBl0xJrw6uvErUvigm/TXBvYE7g4wO1asH331u9IyEhEBrahBEjRmRZDJlJtK8CC4DmqjpLRMoBeXFuon2cm90G2M+XUKimFzDfiecyDMOuc6XO3Ii/wb+2jUHr1kU//JDR7z1K7RK13R2aw2bPvrm8jSocOADnz8+gdu2su+mX7nLjtzUQqQWUBrap6k6nBSISABwCyqnqWfu2QkAkEAwcBLqmVarR1KM1jDv34o8v8t6693iyTAcmDliIT+Wq5F+3KdvPGgsJsZLrrcqUgYwuuuDqpWySnqioiKwC1gPzgO0islBE8tzpyZNS1ThVLZSQZO3bTtrrKNxjf/bYouOGkd293fJtiucuzucHFrD83y3Iv34LfJPh2yUe62CKVVjiOXDgZ/74448siSEzXQfvAzuwlrTJB+S3bwt3fliGYWS1VQdXceXGFbxt3jxbZA2HnnoUqlVzd1gOC05xAWDBZnuCqVOnZkkMmUm051S1t6r+qarn7Y/FwHlXBWcYRtaI2hdFt7nd+Lbbt8QMieG7HguoWXklUb5H3B2aw8aOBXuFxEQBAcLo0T97ZKK97QJcRHyAus4LxzAMd1h/ZD2RXSIJKxtG4YDChJUN4+V6L0FxV8UAACAASURBVLP+9++s6bnR0e4O8Y717AkzZlh9sgneeQeGD69OnjxO6flMV2YSbZy96HdzEWktIv8HbAZMsW/DyOaGNhhKWNmwxPfHLx4nfEU4O49tQ7/9Fl7J3rPfe/a0bnz9/rv1PnduuHTpElOnTuW3335z+fkzk2jfxpqo8BXwP2Ak8D3wsgviMgzDjYrmLsrYZmP56PRSNv/rIWuM1OrV7g7LYbVqWZMXfvgBfHx8GD58OAsXLnT5ee9keJcARYBYzWxjFzPDuwzDeeI1nqafNmX3oU0cnB6Ad4lSsG6dNeo/G5swAfz94YUXIDY2lozMKHXlcuO5UtquluMJSVZE+tzpyQ3D8Fw2sfFJx084632dSR0Kw4YN8MUX7g7LYUOHWkkWyFCSdYa0fjXtEJHeCW9EJEZEbtz6ACJcH6ZhGO5QvmB5pj48lQLPPI+OG2ctxpUDnDsHW7bAyZMnGTBgACtXrnTp+dJaYeFxkt/o+hD4GzgKJHQZ2ICurgnNMAxP8EyNZ6wXtUFVyd7zxCzdu8PevbBxY27mzJnDAw88QOPGjV12vlQTraquuWXTj8A/qnoy6UYR2eiKwAzD8CyzNs9izU8f8d/515DPPoPy5d0d0h17+GEYOBAOHfIjNjYWm4v7nTNz9O+A2/5uUNXTzgvHMAxPdfXGVb47upxrf66Hl7P3YKOEHpBFi3B5koXMJdp3sLoOkhGRzs4LxzAMT/V09aepUasNYxrGw3ffwdKl7g7pjpUrB5UrW8O89uzZQ8uWLV3aT5uZRNsI+E5EokRkmf2xAvgkvYaGYWR/IsLM9jOZ2TiQw4VzoS++CNevuzusO9auHaxYAT4+hTlx4gQXLlxw2bkys9z4MawyhklvhgH4ODUiwzA8Vok8JZjY4b+8sOUJ5kVuhVmzrLVisqH+/eFf/4LSpfPx119/ufRcmUm0k7Fuhl1L2GBf9eBbp0dlGIbH6nF/Dyq+cy903AVdurg7nDuWtPZBAlVFXFB/N8NdB6q6BQgUkdIiEiwiwUBxnDAFV0Tyi8hcEflHRLaLSD0RKSgiS0Rkl/25gKPnMQzDcSJCaMna0KMHW07v4NqVS+4O6Y79/rt1Zbtmze8EBwezbt06l5wnM4W/3wROAvuBffbHBqCJE+J4D/hRVe8DqmGN330FWKqq9wBL7e8Nw/AQO07soM9r1Thbughs25a4kGN2Wmts716IiICTJ4OpW7cuvr6+LjlPZroOamKVRGyEtSLuaaAFsM2RAEQkL9AYeApAVa8CV0WkI9DUvttnwHLgP46cyzAM56lYuCJB9zfA64tVRPfuSqm124jq1YRuc7sR2SXS3eFlSKtW4OUFa9YUJzLSdTFnZtTBSlX9HfgIaKyq+4GPAUd/dZUDYoFPROQvEflQRHIDxVQ1BsD+XNTB8xiG4WQf/ns+E1r4U2qtdb2VkGSTllz0ZAUKQMOG1jAvgPPnzxMfH+/082Qm0VYUkZlAAaylx98A/gtUdjAGb6yr5QhVrQFcJBPdBCLSR0Q2iMiG2NhYB0MxDCMzCo6fwls/3OyjjR16grByzRK7EbKDdu1g82b46KMFFChQgG3bHPojPUWZSbRDsLoLFBiDtRJuA+BFB2OIBqJVNaEXei5W4j1mH9WQMLrheEqNVXWGqoaqamhWVeIxDMMuPJxle5bycisrlTQfmJ+ovcuyXaINCYF8+WowbNgw8ubN6/RzZLoe7W0HEKmgqrsdPMavwL9VdYeIhAO57R+dVNVxIvIKUFBVh6Z1HFOP1jCyVsJaY3NbfUKTau2J2rOUbt8+lq26DxJSYFqjuhytR5vhm2EiEgg8BgRx80pYsG6INbrTAOxeAGaLiC+wF3jafo5Ie6nGg5gqYYbhcRLWGmtSNgx9/XXy+OUlsksk64+szzaJNiHB3rgB16/fYOfObVStWtWp42kzM+pgCVAea8nxpPPuHC7ho6obsZYxv1VzR49tGIbrDG1w84/MjzqW5mDP2gy4tydh02e5MarM27YNmjaFzp1n8sEH/di7dy9ly5Z12vEz00dbBghR1UaqGpbwALo7LRrDMLKtLpW7UOKyNwU/ngOHD7s7nEypUAEuXYKzZ9vwxRdfULBgQacePzOJNgKr2+BW5la/YRjk98vPseeegPh4zo8b7e5wMsXX1xpTu2JFMD17PkG+fPmcevzMJFoBpovI60keI4EPnBqRYRjZ1r86vsbXVcH3o0/gdPYqVd2uHRw5AkuXHnP6yriZSbQNgGJAMyDM/mgOVHFqRIZhZFvlCpTjzydakOvSVXTaNHeHkylt2lg3xt5++3M6dOjAsWPHnHbszNwMm4hVeyDZtAkRcXTEgWEYOcjQfrO4fm0m3g9nr4Ucixa1liKvWLEHo0c3oVChQk47dqbG0YpIc6ypsV+KSEPgqKNjaJ3JjKM1DM9xPf46NrFhE9cvFeNqjo6jzUz1rpexhnj1BFDVVcCLIvLwnZ7cMIycadfJXbR4NZi9T3fKVqswqMLq1TB79kamT5/utONm5ldNF6y6Bj8n2fY18K7TojEMI0coW6AsFY5fp8JnC8CFVbFc4bHHYMyY7xkwYIDTlrfJTKL9n6r+c8u2GoDzOjIMw8gRvG3eVHn6P/xdBOLGhN+c5+rhRKzRBwcPDuDQoeMEBgY65biZSbRxInIP9vXCRKQLMBqY55RIDMPIUXqHPsvUJn4EbN8FP/7o7nAyrF07iIsrzObNzlvUJTOJdgrQHxgmIpeAWUAkMNhp0RiGkWPkzZWXvL2e41BeuPLmKHeHk2HNmoG/P0yZMp9Ro5wTd2bWDLuiqoNVNQgIBgJU9VlVveiUSAzDyHFeaDSEMy/1x7dGbbh2Lf0GHsDfH5o3h3XrfuXjjz/mxo0bDh8zw8O7RKQMVkHu/6jqORFpCpRU1dkOR+EkZniXYRjOcOQI5M59mbx5cyEiWTe8C2vdrtrYJzmo6nKgkIgMv9OTG4aR812Pv86wJa/wvxlDYc8ed4eTISVKQL58fk4rlZiZRLvVvpLBqSTb1gADnRKJYRg5krfNm7+2LyNswETix4xxdzgZFhEBDz74Di++6OgiMplLtClV6foX1hpfhmEYqerT8hVm1lCYNQsOHXJ3OBkSEwPr1x9m+3bHr8Izk2hXiMh3IjJARIaIyHJgABDucBSAiOwXkS0islFENti3FRSRJSKyy/7svPEWhmFkmY4VOxLZuhTx8ddh8mR3h5Mh7doBTOLJJx2v5JWZUQfLgRHAfVgVvHYADVT1C4ejuClMVasn6XR+BauQzT3AUjKxOq5hGJ7Dy+ZFtzYv82VVuDH9Azh1Kv1GbhYaahWa+eabGICKjhwrU9UeVHWbqg5Q1baq+pyqrnHk5BnQEesmHPbnR1x8PsMwXOTp6k/zzzPt0dwB1toxHs5mg7ZtYeHC5gAOTRFLc3iXiHRI8nabqu4WkUrAp1gZfh7wvKpeSql9pgIR2cfN5cynq+oMETmjqvmT7HNaVQvc0q4P0AcgODi41oEDBxwNxTAMV7p2DXx83B1Fuvz9/bl8+XLie1W94yEI6V3Rfg+8CMQAe0TED/gBq7jMeOAyMO5OT36LBqpaE3gYeF5EGmekkarOsI+GCC1SpIiTQjEMw1W2ndnF3E1fwW6PqbCaorff3ouXVw8gwOFjpVf4+xjQTlXjAERkEFAWaKuqi+3bvnQ4CkBVj9ifj4vId8CDwDERCVLVGBEJAo4741yGYbjP2F/H0umNSOLPlMC2e7fHXt1OnBjEjRt5sa4nHRtPm94V7U9JkmwB4D/AooQka+dwHTERyS0ieRJeAw8BW4EFQC/7br2A+Y6eyzAM9xpcdzCfVrmO7eBBqyahhzp4EKxrzb5AJYeOlV6i9U3yegKQG/i/W/ap61AElmLAKhHZBPyOlcx/xOqWaCkiu4CWOK+bwjAMNwktEcqFFo3YUdwHvvvOY0soBgeDdRtqGuDv0LHSS7S/icgCEZkN9AZeU9UdACLiKyITcMLijKq6V1Wr2R9VVHWsfftJVW2uqvfYnz1/TIhhGGma8NsEWpR/iDH17EVm5swhal8UE36b4N7AbjF2LAQ43j0LpNNHq6rTRCQGa9xsT1WdAyAi5YFuwA2sK13DMIwM6Rr5N2WnfH5zQ8+ehAEhA5+01tr2ED17Ws8jRoCjg5kytTijpzPVuwwje4jaF0W3ud2IHXqC4m8VYs5j3xBWNszdYaUqK6t3GYZhOEVY2TD6hfYDoGmFFoTmKgtbt7o5KtdJNdGKSC4Rcbj/1TAM41ZR+6KI2BDBosdD+frvrznctCZ07gxJJgjkJGld0X4BfJVVgRiGcXdI6DaI7BJJ2y/X8+h9jzKg7mnYuRPefNPd4blEWok2CGiU8EZEUqwzYJ9IYBiGkSHrj6wnsktkYp/s550+5+8HivNVzVzouHHw999ujtD50kq0S4HzSd4/mMp+rZ0XjmEYOd3QBkOT3fjK7ZubWY/OYkDzK8T5e0OfPhAf78YInS+t4V2RwCYRCQTigQIi0v2WfXywrnw/cVF8hmHcBZqXa07vlkPZlO8c9dfHwIULkDevu8NymlQTrapuE5GaQEOgJPAY8M0tu/kAT7guPMMw7hbjW46HFgpOWqfLk6Q3YeEqsAxARG6o6m0FZETkdxfFZhjGXUaB/66biu/eAzy3Jz+8+qq7Q3KK9Kp3JUpIsiJSGCgF7FHV86q62VXBGYZxdxER1h1eR+lP5sCSeKhaFR7J/vX+MzxhQUQC7DUPjgF/YJUwnCoiGU7WhmEY6Xm39bt81qwgu0r6owMGwLlz7g7JYZlJklOwqmx1BXbZtzXEWpzRY6/vr127RnR0dLJK6YZr+Pn5UapUKXw8tL6okT0UDijMe+3f54k93Vj70REYPhz++193h+WQzCTaPKra4pZtW0Qk3InxOF10dDR58uQhJCQEyYGd7J5CVTl58iTR0dGULVvW3eEY2VyXyl2Y06ITEVsX0O/995EnnoC6zqjI6h6ZSbTbU9lezRmBuMrly5dNks0CIkKhQoWIjY11dyhGDiAivN/2fXZW7o1ErrX6arOxzCRam4g8DazEGtZVFegHHHJGICLiBWwADqtqOxEpizUFuCDwJ/Av+yiIOzm2M0I00mF+zoYzFQ8sTvH728L9bTl16RQFNfsO/cpM9a6xQD2sK9u/ga+x1vB6wUmxDCL5VfN4YLKq3oO1Om5vJ53HMIxs5NONn9L+lTJcqX4/7NqVfgMPlOFEq6pXVLUP1kywekAJVX1cVc+n0zRdIlIKaAt8aH8vWMXG59p3+QzIkjEes2dDSIi1pntIiPXeUV5eXlSvXp2qVavSvn17zpw5k+zz1q1bc/jwYQBiY2Px8fFh+vTpyfYJCQnh/vvvp3r16lSvXp2BAwcmfnb9+nUKFy7MsGHDkrVp2rQpFStWpFq1ajRo0IAdO3Y4/mUMI4s1K9uMYwHxXNu9A+3b12OXvkmTqrr9gZVQawFNsZYzLwzsTvJ5aWBrKm37YHU5bAgODtZbbdu27bZtqZk1SzUgQNX6l7QeAQHWdkfkzp078fWTTz6pY8aMSXwfFxentWvXTnw/bdo0bdiwoTZp0iTZMcqUKaOxsbEpHn/RokVav359LVeunMbHxydub9Kkia5fv15VVadPn67t27d37ItkQGZ+3oaRUe///r72aWf/n/LTT7P8/MAGdSDHub3wt4i0A46r6h9JN6ewa4q/xlR1hqqGqmpokSJF0j1f06a3P95/3/ps2DCIi0u+f1wcDBpkvT5x4va2mVWvXr3Eq1eA5cuX0zTJgebMmcOkSZOIjo5Otl9a5syZw6BBgwgODmbt2rUp7tO4cWN2796d+YANwwM8F/ocOzs1Zk0ZL+IHvwgvv+zukDLF7YkWa5WgDiKyH+vmVzPgXSB/kskQpYAjrg4kOjrl7SdPOuf4N27cYOnSpXTo0CFx2+LFi2nd2iqAdujQIY4ePcqDDz5It27d+Prrr5O1DwsLS+w6mDx5MgCXLl1i6dKltGvXjscff5w5c+akeO6FCxdy//33O+eLGEYWs4mN2qXr8MIjueD8BZg4EcAjF3VMUUYvfYGyKb125gN714H99TdAd/vrD4D+6bWvVavWbZf8mflTtkyZ5N0GCY8yZTJ8iBTZbDatVq2a5suXT5s1a6bXr19P/KxmzZp65coVVVWdMGGCDh8+XFVVN23apKGhoUliS7nrIDIyUnv06KGqqidOnNBSpUolHr9JkyZ67733arVq1bRjx4568OBBx75IBpiuA8NVlu1dpoXGF9INc6eqgi7bu0wLTyisy/Yuc/m5cbDrIL3ENxHojDUj7P0k29935KRpnC9poi0H/A7stifdXOm1dzTRurqP9syZM9qwYUN97733VFV1z5492rFjx8T9atSooSVKlNAyZcpomTJl1MfHR3fu3KmqqSfaTp06adGiRRPb+Pv765IlS1Q1eR9tVjGJ1nCZkSNTvhIaOdLlp3Y00abXdfAGcAFrCFdTEZkhIq8DwfbiMk6lqstVtZ399V5VfVBVK6hqV1W94uzz3apnT5gxA8qUsYbrlSljvU9YdthR+fLlY8qUKUycOJFr164l6zbYsWMHFy9e5PDhw+zfv5/9+/czbNgwvvoq9dWEzp07x6pVqzh48GBim2nTpqXafWAY2Vp4OKgy8H/WiNKjBXNBbKy13cOlmWjVqs71k6q+CqxQa3jX51gTHQaJyHQRmZYVgWaVnj1h/36rwPv+/c5Lsglq1KhBtWrV+Oqrr/jxxx8TE+2cOXPo1KlTsn07d+6cLGkm7aN98sknmTdvHs2aNSNXrlyJ+3Ts2JEFCxZw5YrLfy8ZRpaL2hfFnK3W/xPFT11h67ghbo4og9K63AUGA3WwZoIl7TqIcOQy2lUPR7sOstLly5c1pXizO0/9eRvZX9I+2eP/11+jyqD78osu27XE5efGxV0Hi4HKwDtYXQfvisgzQCH7lFnjDuXKlYsNGza4OwzDyDaSLupY5O1pbOxcn5AzyuHID90dWrrS6zr4R1U/UdUXgBXAcKxpssWBKSIyU0TezoI4DcO4y926qGOrIREcy2uj7ZIDbowqYzJVvUtV44A1IjJXVacAiEge14RmGIaRukolHuDGh3PwKlbc3aGkKzO1Dqak8trhWgeGYRh3wqtrN642qMu66HXuDiVNnjAzzDAM44699Xkf1nVvyPHY/e4OJVUm0RqGka09U6gFA3+7TtS459wdSqpMok1iwm8TiNoXlWybM+ZSe3KZxP379+Pv70/16tWpVq0a9evXT9xv+fLl5MuXL/GcLVrcupKRYbhf6fY9iS6dj3u/WsLR8zHuDidFJtEmUbtEbbrN7ZaYbKP2RdFtbjdql6jt0HH9/f3ZuHEjW7dupWDBgkybdnOOx6VLlzh16hQlS5YE4JtvvqFu3bopzu6Kiopi48aNbNy4kSlTErvJ+fnnn6lYsSKRkZEJ45wTzZ49m02bNtGrVy9eTqXiUfny5dm4cWPifm+++WbiZ40aNUo85y+//OLQz8EwXEIEnxdepMYR5euPXnJ3NCm665YKb/pp09u2davSjf61+1OnVB1K5ClBq1mtCMoTRMz5GCoVqcSBs9bwkRNxJ+gS2SVZ2+VPLc/U+evVq8fmzZtvtk+lTGKPHj04fPhwYgJOS0KZxIiICNauXUu9evVu26dx48a8++676R7r3LlzFChQIGNfxjA8RLG+Q7g08k0qfr2U+EHx2MSzriHvukSbngJ+BQjKE8TBswcJzhdMAT/nJZ2EMom9e99clWfx4sU88oi1eERKZRJfeunmb+iwsDC8vKx5Ir169WLw4MGJZRKnT5/OmTNnmDNnToqJNq0yiXv27KF69eqcP3+euLg41q27eQf3119/pXr16gB07dqVESNGOP6DMAxny5MHr/7P0+rSVSTFctZu5si0Mk97OGMKbsI0v9eWvea0EmyeXCZx3759WqVKlcT3X331lbZq1UpVVaOiorRt27aZ+q5mCq7hbqfiTumJiyecekyy+woLniShTzaySySjwkYR2SUyWZ/tnUrooz1w4ABXr15N7KPdu3cvpUuXxtfXF7C6AD799FNCQkLo0KEDmzZtYlc6i9HNmTOHX375hZCQEGrVqsXJkyeJiroZ7+zZs9m4cSPff/89pUuXZt26dYk3txYsWHDb8Tp06MDKlSsd+r6G4S4Xrpyn29CyjF4W7u5QkjGJNomkc6kBwsqGEdklkvVH1jvl+J5QJrFOnTqJN7eSrvSQYNWqVZQvX97xL2sYbhC4Yg1L/nuWQ3M+4PC5jC0FlRU8oo9WRPyAlUAurJjmqupIESmLtbxNQeBP4F+qetVVcQxtMPS2bWFlw5LNr3bUrWUSp06dCqReJrF79+689tprVixJ+mgfeOABmjVrlmKZxKFDh2aqTGJCH62q4uvry4cfen6RDsNIUbNmXA8qznPrjjH+t/FMeXhK+m2ygiP9Ds56YC3GGGh/7QOsA+oCkSRfzqZfWscxZRLdz1N/3sZdZNQoVdAqg3w0+my0Uw5JTuijtX+XC/a3PvaHYi3UONe+/TPgETeE5xKmTKJhuMizz6I+Pjz7+w0W7lzo7mgAD+qjFREvEdkIHAeWAHuAM6p63b5LNHDboFIR6SMiG0RkQ2xsbNYFbBiGZypeHOncmeePlKRvzT7ujgbwoESrqjdUtTrW0uIPApVS2i2FdjNUNVRVQ4sUKeLqMA3DyA7eeQfvbf+AzcbpS6fdHY3nJNoEqnoGWI7VR5tfRBJu2JUCjrgrLsMwspGgIAgI4PONn1HqnZIcOnvIreF4RKIVkSIikt/+2h9ogbWSQxSQMOe1FzDfPREahpHtbN3K411HUWv/Vd789c3093chj0i0QBAQJSKbgfXAElX9AfgP8JKI7AYKAR+5MUbDMLKTkBB8jp9g4q6yfPTXRxw4474lbzwi0arqZlWtoaoPqGpVVR1l375XVR9U1Qqq2lVVs24NbSeuFe/JZRI7derE999/n/i+YsWKjBkzJvF9586dmTdvnimZaGQ/gYHw1FPUXn2AYhdx61WtRyRaj/TGG047lCeXSaxfvz6rV68G4OTJkwQGBrJmzZrEz9esWUP9+vUBUzLRyIb690euXeOtfeWZvWU2Zy+fBZxTZzoz7r5E27Tp7Y/337c+i4u7uS3pvp9+ar0/ceL2tplUr169xKtXSL1MYnR0dLL90pJQJjE4OJi1a9emuE/jxo3ZvXv3bdsbNGiQmGhXr15Nu3btiI2NRVXZt28f/v7+FC/u+YvfGUaKKlbkVINahP20k4/bzCCfXz6n1ZnOjLsv0aZl7FhYscJ6wM3XSf60dkRCmcSkNQaS1jtIqUxiUmFhYYl/uk+ePBkgsUxiu3btePzxx1Otc5BamcRatWqxdetWrl69yurVq6lXrx4VK1Zk+/btrF69mgYNGiTum1AysXr16owdO9bhn4dhZIWCb03m9NhXGfjjQF5b9hpdv+marKZJlnBkWpmnPZw6BTdxwprjPLlMoqpq/fr1dc2aNdq0aVM9deqUTps2TWfOnKn9+/fXiIgIVc14yUQzBdfwVK8te00JRwtPKKwXr17MVFtywhTcnM7TyyTWr1+flStXcv78eQoUKEDdunVZvXr1bVe0hpFd/bpxAf7jJvHTuns5EXeCx+Y+lqXnN4k2NSNHOv2QnlomsUGDBkyfPp1q1aoBVmWwtWvXcvDgQapUqeLEn4BhZL2ofVH0+e5p/hN1lYcW7+SxKo/xw84feGO58254p8ck2tQ4cXhXUreWSUxItKmVSUyaNJP20T755JPMmzcvxTKJCxYsyFSZxPr167N3797EJXC8vb0pWrQooaGh2GzmPxEje1t/ZD3vPzUXW9duAHxR+VXuK3wf434bx8GzB7MkBlG9rXxAthUaGqq3VsTavn07lSqlVDbBva5cuUKDBg1yXAUvT/15G3ex8PAUh2tOfbggodN/oF7p29fYu5WI/KGqoXcagrlccRNTJtEwskh4OFi3t6339epBhQr0n3swQ0nWGUyiNQzj7rJiBfzyC14Bubked4FJX7/IygOuXSfPJFrDMO4eI0eCjw+UKQNA/H/+wzPPTGXm2Ec5demUy05rEq1hGHePW25y+w4aTK6QCnwx8yRLe9ZHr19PuZ2DTKI1DOPuVaECARs2sqX9g3T9bgeH61WFY8ecfhqTaA3DuLv5+1Nl/hreebYqvtt3cu6484uEm0SbgpiYGJo0acLRo0edcrzAwEAA9u/fT9WqVW/7fO3atdSpU4fq1atTqVIlwl00htcwjJTZxEbPd34hetNK8t4fao1QWLDg5kgFR4/vlKPkMKNHj2bVqlWMGjUqS87Xq1cvZsyYkVhKsVu3bllyXsMwbioWWIya5RsCsPfLadCxI3TqBGfOUApKOHJstydaESktIlEisl1E/haRQfbtBUVkiYjssj8XcMb5mjZtyqf2sofXrl2jadOmzJo1C7BqEogIERERxMfHExERgYgk1iI4ceIETZs2ZeFCawljZ13xHj9+nKCgIMAqEl65cmWnHNcwjMz7bvt3lN/1ApuGPQ2LFkHNmhSzVoG5Y25PtMB1YIiqVsJakPF5EakMvAIsVdV7gKX29y61detWihYtmphY/f39KVq0KBERES497+DBg6lYsSKdOnVi+vTpXL582aXnMwwjdW3vbUtoyVCa5vuOo/+LhGvXHD+oI6W/XPHAWoCxJbADCLJvCwJ2pNfWGWUS+/btqzabTf38/NRms2m/fv0y1T4luXPnVlXVffv2aZUqVVLcZ/fu3fr+++9r48aNtUmTJg6f011MmUQjJ9h1cpeObeaTMJ9Ma1llU3NGmUQRCQFqAOuAYqoaA2B/LppKmz4iskFENsTGxjocw7Fjx+jbty9r166lb9++TuseSE/58uXp168fS5cuZdOmgaQrmgAACyFJREFUTZw8eTJLzmsYxu0qFKzA5n6PIuHwRlS4w8fzdvgITiIigcC3wIuqek5EMtROVWcAM8AqKuNoHPPmzUt8nXRtL1datGgRbdq0QUTYtWsXXl5e5M+fP0vObRhGyp6r9Rzzd8zn0rVLDh/LIxKtiPhgJdnZqpqQ6Y6JSJCqxohIEHDcfRE6z44dOyhVqlTi+8mTJ/Ptt98yePBgAgIC8Pb2Zvbs2Xh5ebkxSsMwwsqG8b8e/6Pb3G7k8uWCI8dye6IV69L1I2C7qr6T5KMFQC9gnP15vhvCc4oLF6x/o5CQEK6l0LHetWvXrA7JMIwMCCsbRr/QfowuNDrQkeO4PdECDYB/AVtEZKN923CsBBspIr2Bg4DJRoZhZKmofVFEbIiAi8Q4chy3J1pVXQWk1iHbPCtjMQzDSJCwLHlkl0iaDW12xJFjedSoA1fRHLSKhCczP2cjJ1l/ZL3TliV3+xWtq/n5+XHy5EkKFSpERkcyGJmnqpw8eRI/Pz93h2IYTjG0wVCnHSvHJ9pSpUoRHR2NM8bYGmnz8/NLNqLCMAxLjk+0Pj4+lC1b1t1hGIZxF7sr+mgNwzDcySRawzAMFzOJ1jAMw8UkJw3JEZHzWFW/cqrCwAl3B+FC5vtlb//f3rkHX1VVcfzzBQTBHyIQIcIvQQKMwUCHIV429EBBHbFJIiRBBzOJEpiaRhJTi8JHKToiFaEUITAaJlGjEwg5iUNCJmKoaGBCvMxAcpzisfpj7xvHy72/1/1dzu+eWZ+ZPffufc7eZ6277qzZZ52z186yfn3MrG1DO2ftYdirZjYwbSHKhaSNrl/l4vpVLpI2ltLfQweO4zhlxh2t4zhOmcmao/1p2gKUGdevsnH9KpeSdMvUwzDHcZymSNZmtI7jOE0Od7SO4zhlJjOOVtIoSa9Kel1S2bcmLzeSqiWtlbRV0suSpsX2DpJ+L2lb/GyftqwNRVJzSS9IWhXrPSRtiLotl9QybRkbiqQzJD0m6ZVowyEZs92M+L/cImmppFMr2X6SHpK0T9KWRFtBeylwf/Q1myVdUNv4mXC0kpoD84DRQF9gvKS+6UpVMkeAb5jZx4DBwNSo003AGjPrBayJ9UplGrA1Ub8TuDfq9i9gcipSNQ73AU+a2blAf4KembCdpK7AjcBAM+sHNAe+SGXbbxEwKq+tmL1GA71iuR6YX+vopexV3lQKMAR4KlGfCcxMW65G1vEJYCRh5VuX2NaFsEgjdfkaoE+3+Of9NLCKsMvG20CLQjatpAKcDmwnPmxOtGfFdl2Bt4AOhEVPq4CLK91+QHdgS232An4CjC90XrGSiRktxw2fY2dsywSSugPnAxuAzma2GyB+fjg9yUpiLvAt4FisdwQOmNmRWK9kG54D7AcejqGRn0k6jYzYzsx2AT8k7OW3GzgIbCI79stRzF719jdZcbSFtk7IxHtrkqoIW7FPN7N305anMZB0GbDPzDYlmwucWqk2bAFcAMw3s/OB96jQMEEhYqxyDNADOAs4jXA7nU+l2q826v1fzYqj3QlUJ+rdgJI2U2sKSDqF4GSXmNmK2LxXUpd4vAuwLy35SmAYcLmkHcAyQvhgLnCGpFz+jUq24U5gp5ltiPXHCI43C7YD+Cyw3cz2m9lhYAUwlOzYL0cxe9Xb32TF0T4P9IpPPVsSAvMrU5apJBQ2OFsIbDWzexKHVgKT4vdJhNhtRWFmM82sm5l1J9jqaTObAKwFroynVaRuAGa2B3hLUp/Y9Bngr2TAdpG/A4MltYn/05x+mbBfgmL2WglMjG8fDAYO5kIMRUk7AN2IgexLgNeAN4Cb05anEfQZTrgd2Qz8JZZLCLHMNcC2+NkhbVlL1HMEsCp+Pwf4E/A68CjQKm35StBrALAx2u/XQPss2Q64HXgF2AIsBlpVsv2ApYR482HCjHVyMXsRQgfzoq95ifD2RY3j+xJcx3GcMpOV0IHjOE6TxR2t4zhOmXFH6ziOU2bc0TqO45QZd7SO4zhlxh2tkwoxc9fkmP1oRAP6t5b0pqShZRCvbBTTOy7VHZuiaE4ZydouuM5JIjq47xJWdf0SeB8YBDwJ3GLH17wXxMyOStoAnNeQ65vZ+5Jm88HsX42GpCuBscC7hOxUbYE7zayk3VBr0PsuQi4LJ4O4o3UahJmtl7QcGGFmEwEk9Se8pH8MuLkOw7xcn2tK+raZ/SAhw4L69K/Hde4mOMKxZnYotvUCVkuaYceXQzeUE/Q2s6V5MgwBWprZH0q8ltME8NCBUwqHkxUze5GwUujyunS2eqyWkXQLcFG9pGsAksYAM4DJOScLYGbbgDnAYklnlXKN2vSWdDawhMLJS5wKxB2t09i0Jy9ZiqRhkmZJ+oWkdZJ6FuooqVM8Z5akpyTNie1DCbl4e0iaLek8SV+V9JKkEZI6SnpEkkm6IfY5XdKyGF5AUjdJt0q6R9KfJV1cRP5pwCYLqQDzeRxoA1wnqb+k5ySti+P3lfRMot5M0o8kfT/q9JuYKrGQ3hMkrZc0KSZlGU/IjDVZ0tSosynsAtAm9rle0kZJ1YXGdJoYaa8x9lK5BbgGOBK/NyOEC94DPpk4pxuwIFFfDjybqBsh/ABwN7Awfv9EPNYx1m8D1iX6Vef1rQIOAJ9LnHM/YVbYjJAFrVVsn0KIvXYsoNNBYFkNOh8Cnigi0//rwKXAG/F7i9jv80X0rorXvabQ8VhfAzyYqF8HjEr7P+ClbsVjtE6pNJN0O/AlYC/Qz8y2J45PANpJmh7r+4EqSc3N7GjeWPMJicvaAxfGtirgnwWuuzNZMbN/S1oEfAV4XNJHgB1mZpIGEWaIU0KyKToA6wnJmvPHbgnky5XkICGBSm08A0yI2eQuJYRZqgqdGGU/UMt4s4FVkr5jZm8TZvjj6yCH0wRwR+uUiggz0b8BCwiJoJOO9mxgs5nNrcNYu4FZhC1Rnk6MfwLRgeY3PwhslfRRYBxhy5GcDO/UUYYdwJmFDsSUgO0J2alqxMwOSfo4IUH2QsIMuqaYa41xWzNbK+lFYLqkXwEvmNmxmvo4TQeP0TqNgpn9HHgIeFRS0lH9A/iCwgaaAEgargJekpB67piZ3Qu80wAZXiM46K8DbePMLyfDsDjLzcnQWVLvAsOsAAblYqF5DCDEaJfE+lHCq18nIGkicLWF3Lu1OuY6MgeYSogjP9xIYzonAXe0Tink7ohyzuZGwl5KKySdGtuWAb0Jt72jJI0DRscZac7Z5j4HAp2iUx4Z2zrGp/z/ATrEF/77FOibYx7wNcKGgTk2EGbZqyWNkzQSuIMPzrxz3EFwzN9MNsbrzQLuM7PnYvMe4FxJH5LUj7BzRKcYLhgItFXYhnso0A5oI6lnEdmVV/9v1D25m/MqQtLt1ma2t4DsTlMl7SCxl8oshBjqWsIt761AdWyvJrx1sBG4NrZdREiQfJCwuKEqtk+J/RcAnQlZ7A8Q4qfDCU5lMXAKYRv5XcBvCbfvN8S+80kk0CY4/d8VkLcnsJrwUOqPQO8adDuT4NQeIOwecAXhId6svPNaA89GfW+Lv8Nq4DLC7PdNQnLoqwgz5ecJDwfz9R5DcKwrgK5x7AcIDv+qvGveBHwqbft7qV/xxN+OUwRJwwhvK3QGJpjZIymLhKTFZnZ12nI49cMdrePUQIw3/xgYQnjdqixLfmuRYQBhht8O2GVmi062DE5peIzWcWrAzPaY2RWE8MGXJV0rqUH5GUrgQuB7hNjsopN8bacR8Bmt4zhOmfEZreM4TplxR+s4jlNm3NE6juOUGXe0juM4ZcYdreM4Tpn5H56tapykuRpWAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import csv\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import os\n",
    "import re\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "\n",
    "DIR_NAME = './flow_size/'\n",
    "FIG_SIZE = (5, 4)\n",
    "fontdict = {'family': 'Times New Roman', 'weight': 'normal', 'size':15}\n",
    "\n",
    "filenames = os.listdir(DIR_NAME)\n",
    "filenames_dict = {\n",
    "    'R/AEAP': [],\n",
    "    'R/AEAP-BF': [],\n",
    "    'R/AEAP-WF': [],\n",
    "    'D/AEAP': [],\n",
    "    'ILS': []\n",
    "}\n",
    "\n",
    "for filename in filenames:\n",
    "    if re.match(r'b.*(ba|rrg|er).*backtracking.*aeap\\.csv', filename):\n",
    "        filenames_dict['R/AEAP'].append(filename)\n",
    "    elif re.match(r'b.*(ba|rrg|er).*backtracking.*aeapbf\\.csv', filename):\n",
    "        filenames_dict['R/AEAP-BF'].append(filename)\n",
    "    elif re.match(r'b.*(ba|rrg|er).*backtracking.*aeapwf\\.csv', filename):\n",
    "        filenames_dict['R/AEAP-WF'].append(filename)\n",
    "    elif re.match(r'b.*(ba|rrg|er).*dijkstra.*aeap\\.csv', filename):\n",
    "        filenames_dict['D/AEAP'].append(filename)\n",
    "    elif re.match(r'o.*(ba|rrg|er).*backtracking.*aeap.*\\.csv', filename):  # only care about R/AEAP\n",
    "        filenames_dict['ILS'].append(filename)\n",
    "\n",
    "objective_value_dict = {}\n",
    "DAEAP_objective_value_list_list = []\n",
    "num_scenarios = 19 * 3 * 6  # all of scenarios\n",
    "no_of_scenarios_mapping = {}\n",
    "for alias in filenames_dict.keys():\n",
    "    objective_value_dict[alias] = []\n",
    "    for i in range(19):\n",
    "        if alias == 'D/AEAP':\n",
    "            DAEAP_objective_value_list_list.append([])\n",
    "        else:\n",
    "            objective_value_dict[alias].append([])\n",
    "for alias, filenames in filenames_dict.items():\n",
    "    for filename in filenames:\n",
    "        with open(os.path.join(DIR_NAME, filename), 'r') as file:\n",
    "            reader = csv.reader(file)\n",
    "            i = 0\n",
    "            for line in reader:\n",
    "                if alias == 'D/AEAP':\n",
    "                    # DAEAP_objective_value_list_list[i].append(float(line[3]) - float(line[7]))\n",
    "                    DAEAP_objective_value_list_list[i].append(float(line[3]))\n",
    "                else:\n",
    "                    # objective_value_dict[alias][i].append(float(line[3]) - float(line[7]))\n",
    "                    objective_value_dict[alias][i].append(float(line[3]))\n",
    "                i += 1\n",
    "            \n",
    "for alias in filenames_dict.keys():\n",
    "    if alias == 'D/AEAP':\n",
    "        continue\n",
    "    objective_value_list_list = objective_value_dict[alias]\n",
    "    Y = [0] * 21  # no. of scenarios\n",
    "    for i in range(19):\n",
    "        objective_value_list = objective_value_list_list[i]\n",
    "        DAEAP_objective_value_list = DAEAP_objective_value_list_list[i]\n",
    "        for j in range(DAEAP_objective_value_list.__len__()):\n",
    "            objective_value = objective_value_list[j]\n",
    "            DAEAP_objective_value = DAEAP_objective_value_list[j]\n",
    "            for k in range(21):\n",
    "                if objective_value >= DAEAP_objective_value * (0.0 + 0.05 * k):\n",
    "                    Y[k] += 1\n",
    "    no_of_scenarios_mapping[alias] = [n / num_scenarios * 100 for n in Y]\n",
    "    \n",
    "markers = {\n",
    "    'R/AEAP': 'o', 'R/AEAP-BF': 'x', 'R/AEAP-WF': '+',\n",
    "    'D/AEAP': 's', 'd-aeapbf': '^', 'd-aeapwf': '<',\n",
    "    'ILS': '*',\n",
    "}\n",
    "colors = {\n",
    "    'R/AEAP': 'blue', 'R/AEAP-BF': 'green', 'R/AEAP-WF': 'red',\n",
    "    'D/AEAP': 'blue', 'd-aeapbf': 'green', 'd-aeapwf': 'red',\n",
    "    'ILS': 'black',\n",
    "}\n",
    "linestyles = {\n",
    "    'R/AEAP': '--', 'R/AEAP-BF': '--', 'R/AEAP-WF': '--',\n",
    "    'D/AEAP': '-', 'd-aeapbf': '-.', 'd-aeapwf': '-.',\n",
    "    'ILS': ':',\n",
    "}\n",
    "tick_spacing = 20\n",
    "x = np.array(np.arange(21) * 5) # relative quality\n",
    "fig_1, ax_1 = plt.subplots(figsize=FIG_SIZE)  # no. of available flows\n",
    "# ax_1.set_title('Relative Quality')\n",
    "for alias in filenames_dict.keys():\n",
    "    if alias == 'D/AEAP':\n",
    "        continue\n",
    "    if len(no_of_scenarios_mapping[alias]) == 0:\n",
    "        continue\n",
    "    ax_1.plot(x, no_of_scenarios_mapping[alias], label='{}'.format(alias), \n",
    "              marker=markers[alias], markersize=6, markerfacecolor=colors[alias], color=colors[alias], linestyle=linestyles[alias])\n",
    "    ax_1.set_xlim(x[0], x[-1])\n",
    "    ax_1.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))\n",
    "    ax_1.set_xlabel('Relative Quality', fontdict=fontdict)\n",
    "    ax_1.set_ylabel('# of Scenarios (%)', fontdict=fontdict)\n",
    "    ax_1.legend()\n",
    "    fig_1.savefig('./fig_rel_qual.png')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}